Critical properties of the band-insulator-to-Mott-insulator transition in the 
strong-coupling limit of the ionic Hubbard model 
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We investigate the neutral-to-ionic insulator-insulator transition in one-dimensional materials by 
treating a strong-coupling effective model based on the ionic Hubbard model using the density- 
matrix renormalization group and finite-size scaling. The effective model, formulated in a spin-one 
representation, contains a single parameter. We carry out an extensive finite-size scaling analysis 
of the relevant gaps and susceptibilities to characterize the two zero-temperature transitions. We 
find that the transition from the ionic band-insulating phase to an intermediate spontaneously 
dimerized phase is Ising, and the transition from the dimerized phase to the Mott-insulating phase 
is Kosterlitz-Thouless, in agreement with the field-theory-based predictions. 

PACS numbers: 71.10.-w, 71.10.Fd, 71.10.Hf, 71.30.+h 



Electrons in solids are subject to both a single-particle 
potential and the Coulomb interaction. A wealth of inter- 
esting phenomena can occur when the form of the single- 
particle potential deviates from that of the ideal crystal 
due to, for example, structural transitions, lattice vibra- 
tions, or defects or impurities. A simple Hamiltonian that 
incorporates the combined effects of interactions and re- 
duced translational symmetry in a particularly transpar- 
ent manner is the ionic Hubbard model (IHM), in which 
the single-particle energy alternates between neighbor- 
ing sites. This model was introduced by Nagaosa and 
Takimotoi'^'^ to describe the neutral-ionic transition ob- 
served by Torrance et al. in mixed-stack organic charge- 
transfer compoundsi^ In a mixed stack of donor (D) and 
acceptor (A) molecules, the neutral phase corresponds to 
a uniform and neutral distribution of charge, D°A"D°A°, 
and the ionic phase to an alternation of positive and 
negative charges, D"'"A~D+A~. The insulating behav- 
ior in the neutral phase originates from the Coulomb in- 
teraction between electrons, i.e., the Mott mechanism, 
whereas the ionic phase is essentially a band insulator. 
Recently, the neutral-ionic transition has been observed 
in organic charge-transfer compounds close to zero tem- 
perature, motivating interest in it as a pure quantum 
phase transitioni^ 

A different class of quasi-one-dimensional materials 
in which a similar charge disproportionation occurs is 
that of the halogen-bridged transition-metal complexes, 
whose structure is formed by a backbone of of alternat- 
ing metal and halogen atomsi^ In these MX-chain com- 
pounds (or in the related MMX materials^), a sponta- 
neous breaking of the translational symmetry occurs due 
to the dimerization of the halogen sublattice, XMX-M~ 
XMX-M. The differing distances of the halogen ions from 
the neighboring metal ions give rise to a two-fold alter- 
nation in the energy of the d levels. 

The Hamiltonian of the ionic Hubbard model can be 
grouped into three terms, a one-dimensional nearest- 
neighbor hopping term with matrix element t, an on-site 
Coulomb repulsion of strength U, and an ionic alternat- 



ing potential of depth A, 

H = Ht+Hu + HA, 

with 
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Here (cio-) are the usual creation (annihilation) oper- 
ators on site i for an electron of spin a and fiicr = cj^ Cio- . 
Without the ionic potential, A = 0, the model reduces to 
the one-dimensional Hubbard model, whose behavior is 
well understood/^ Although the overall physics described 
by the ionic Hubbard model is now fairly well known, 
many details of the transition are still unclear. The gen- 
eral behavior in the ground state is summarized in the 
schematic ground-state phase diagram shown in Fig. [T] 
When A > i7, the system is a band insulator (BI) and 
has both a charge and spin gap. When A < U, the sys- 
tem is a critically antiferromagnetic Mott insulator (MI) 
with a charge gap and gapless spin excitations. These two 
phases are separated by two continuous phase-transition 
lines within which there is a spontaneously dimerized in- 
sulating phase (SDI) of width of order t, i.e., a phase with 
both spin and charge gaps as well as with long-range bond 
dimcr order. 

In order to understand the origin of the phases, let us 
first examine what happens in the atomic limit, t = 0, 
which can be easily treated. For U > A and at half 
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Figure 1: (Color online) Ground-state phase diagram of the 
ionic Hubbard model. Location of the phase boundaries is 
approximate, but drawn to scale according to values from 
Refs. 1^ and [lol . The (red) shaded intermediate region desig- 
nates a spontaneously dimerized insulating phase (SDI). 



filling, there is no double occupancy in the ground state, 
which consists of a series of singly occupied sites with 
energy ±A/2 so that the entire system has energy E ^ 
0. For U < A, double occupancy is favorable, and the 
ground state consists of doubly occupied sites at energy 
U — A alternating with empty sites, so that the energy 
of the system is L{U - A)/2. At (f7 - A) ^ 0, a level 
crossing of two configurations occurs and there is a first- 
order transition at U = A. 

Turning on the hopping term leads to more subtle be- 
havior in the vicinity of the transition. In the nonin- 
teracting limit, U = 0, the Hamiltonian is diagonal in 
momentum space. It follows that the ionic term, A, 
opens a charge and a spin gap, and the two gaps have 
the same value. Correspondingly, spin-spin and charge- 
charge correlations decay exponentially. The scenario 
does not change with the inclusion of a weak interaction 
Hu; the electrons tend to doubly occupy sites with lower 
potential, and the system remains a band insulator. 

In the large-JJ limit, the double occupancy can be 
treated perturbatively and the low-energy physics of the 
IHM is described by an effective spin one-half Heisen- 
berg modelfiiiiii^ It is important to note that this effec- 
tive model restores translational invariance, and that the 
charge and spin sectors are completely separated. The 
system has gapless spin excitations and critical spin-spin 
correlations, while the charge gap, in contrast, scales as 
U for large U. This description is robust for a wide range 
of parameters in the strong coupling limit, but fails close 
to the transition line because perturbation theory breaks 
down in the critical regime i In fact, there are analytical 
and numerical indications that show that higher-order 
spin excitations mix into the charge degrees of freedom 
everywhere in the MI phase 

A few years ago, Fabrizio, Gogolin and Nersesyan 



proposed a new, interesting scenario based on field- 
theoretical arguments They argued that two quantum 
phase transitions occur, an Ising transition between the 
band insulator and an intermediate spontaneously dimer- 
ized phase, followed by (for increasing U / A) a Kosterlitz- 
Thoulcss transition (KT) between the dimerized phase 
and the Mott insulator. This scenario is based on an ar- 
gument in which the transition is approached, on the one 
hand, from the MI limit and, on the other hand, from 
the BI limit. The authors consider the weak-coupling 
case, {U,A) « t, and use standard bosonization. The 
Hamiltonian then consists of three parts, a first term de- 
pending only on charge degrees of freedom, a second term 
involving only spin degrees of freedom and a third term, 
proportional to A, which couples charge and spin degrees 
of freedom. Starting from the MI phase {U dominating) 
with a charge gap but no spin gap, one can integrate 
out the charge degrees of freedom. This leads to a sine- 
Gordon model for the spin degrees of freedom with a 
positive coupling for U > Uc^- The coupling term turns 
negative for U < , and therefore corresponds to a 
KT transition point. A spin gap opens for U < Uc^ and is 
attributed to a spontaneously dimerized insulating phase 
(SDI). Starting from the BI phase (A dominating), which 
exhibits both a charge and a spin gap, Fabrizio, Gogolin 
and Nersesyan calculate spin and bond-order suscepti- 
bilities using perturbation theory. A critical value C/ci is 
found where the bond-order susceptibility diverges, while 
the spin susceptibility remains finite. Thus, Uci must be 
in a region with a finite spin gap, and it follows that 
Uci < Uc2- Close to C/ci it is argued that the spin de- 
grees of freedom can be considered to be frozen. This 
yields a double sine- Gordon Hamiltonian for the charge 
degrees of freedom, which is known to undergo a quantum 
phase transition of an Ising typCf^ The order parameter 
of this transition is the bond order operator, which con- 
firms that the intermediate region, Uci < U < Uc.,, is in 
a SDI phased 

At least one transition has been found in all numeri- 
cal work2^^a2ii8a9i20^^ published after Ref. 0, al- 
though, for the most part, the critical behavior was not 
characterized. The critical exponents were calculated in 
Ref. [13, but were found to deviate from the expected 
two-dimensional-Ising values. However, even confirming 
that there is a second transition has been a quite difficult 
task. The two transitions turn out to be very close to one 
another and, since the transition to the Mott insulator 
is expected to be a KT transition, it is very difficult to 
find and characterize using finite-size-scaling studiesJ^ 
For these reasons, studying an effective model character- 
izing the region of the transition and the intermediate 
phase is useful. 

Another very important subtlety is how to map the 
gaps from the field-theoretical model onto the original 
lattice model. In the ionic Hubbard model, the charge 
gap, the one-particle gap, and the spin gap all behave 
differently at the transitions. The one-particle gap is re- 
lated to the charge and spin gaps but is fundamentally 
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different because it involves a change of the particle num- 
ber, while the charge and spin gaps are spectral gaps of 
excitations into the charge and spin sectors, respectively, 
only. One way of locating critical points is to examine 
the smallest energy gap, i.e., the mass gap, as a function 
of the tuning parameters. The critical point is then the 
point at which the gap vanishes in the thermodynamic 
limit. 

The remainder of this paper is organized as follows. 
In Sec. m an effective spin-one model for the transition 
is derived via a strong-coupling treatment. In Sec. [Til 
the numerical method used to study the model is de- 
scribed. In Sec. mil and IIVI we report the analysis of 
the band-insulator-to-spontancously-dimerized insulator 
and the spontaneously-dimcrizcd-to-Mott insulator tran- 
sitions, respectively. 



Table I: Mapping between the single-site basis states of the 
ionic Hubbard model {|0), | t)) I i)i \d)} with |c?) denoting 
the doubly occupied state, and those of the effective spin-one 
model {jS"^}}. 
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The Hamiltonian for the effective spin-one model can 
thus be expressed in terms of the usual spin-one opera- 
tors, yielding H'^ = + H^, with the exchange term 



I. EFFECTIVE MODEL 

A. Derivation of the effective Hamiltonian 

In order to investigate the critical behavior of the ionic 
Hubbard model at half filling, we derive an effective 
model, formulated in terms of spin-one operators, valid 
for (C/, A) >> t. In this limit, the doubly occupied state 
on the even sites (with on-site potential A/2) and the un- 
occupied state on the odd sites can be projected out. At 
half filling, a double occupancy on an even site is neces- 
sarily associated with a completely unoccupied odd site, 
and has a cost in energy of U + A. This procedure is a 
second-order strong-coupling expansion with parameter 
t/{U, A) analogous to that used to derive the t-J model 
from the Hubbard model. In fact, the resulting model 
can equivalently be formulated in terms of t-J operators 
rather than spin-one operators; we feel that the latter for- 
mulation is more intuitive for the half-filled systemi^i^ 
The physical meaning of the spin-one states is as follows: 
the S'^ = ±1 state corresponds to a singly occupied site 
with a spin-i electron with spin up or down, while the 
Sz = state corresponds to an unoccupied site on the 
even sites and a doubly occupied site on the odd sites. 
The mapping of the states of the ionic Hubbard model 
to those of the effective spin-one model is summarized in 
Table m 

As we shall see, conservation of particle number leads 
to a spin exchange process for the spin-one operators 
that is more restricted than the Heisenberg exchange. 
Given the mapping of states described above, the effec- 
tive Hamiltonian can most easily be derived by first ex- 
pressing the original Hamiltonian in terms of transition 
operators between the fermionic states (Hubbard opera- 
tors), then projecting out the states as outlined above, 
and subsequently writing the Hamiltonian in the reduced 
state space in terms of spin transition operators. Finally, 
the transition operators in spin space can be rewritten in 
terms of spin-one operatorSi ^^'^^ A detailed derivation is 
given in the appendix. □ 



1^1 



(5) 



and the interaction term governed by the single parame- 
ter £ =[/ - A, 



e 2 ^ 



sn -1 



(6) 



Note that it is immediately clear from the effective 
model that the relevant interaction parameter is e = U — 
A. For t = 0, it is clear that there should be a transition 
aXU ^ Is. because the sign of the term changes. For 
£ > > t, the on-site 5^=0 state is strongly suppressed so 
that the remaining degrees of freedom, S'^ = 1 and — 
— 1, correspond to the localized spin-i degrees of freedom 
of the MI phase of the original model. For £ — > — oo, the 
S"^ = ±1 local states are suppressed, leading to a ground 
state that is a simple product of local 5^ = states, 
which maps to the band insulator. However, the nature 
of the transition(s) and possible intermediate phases for 
finite t still needs to be determined. In particular, it 
is important to investigate whether the behavior in the 
vicinity of £ = agrees with previous numerical results 
for the ionic Hubbard model )^°d^'^^ as well as with field- 
thcorctical trcatmentsjii 

Note that the derivation of the effective model can eas- 
ily be extended to include additional interaction terms 
that do not break the symmetries of the original model, 
such as a next-nearest-neighbor Coulomb repulsion. In a 
similar context, a related effective model was developed 
some time ago in Ref. [27l . 



B. Observables 

Since the formulation of the effective model in terms 
of spin-one operators is a notational convenience rather 
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than physical, we are interested in studying observables 
of the original model. Therefore, it is necessary to trans- 
late the observables of the IHM into the language of the 
spin-one model. The local spin operators map as (small 
letters: IHM, capital letters: effective model) 



St - 


2 ' ' 


sf - 








si - 





the local charge operators as 

{^Sf^ i = even 
2 - (^5f j ^ i = odd , 
and total spin and charge operators as 

N ^ L + j2i-^y{s!y ■ (7) 

As we can see, conservation of in the IHM leads to 
conservation of 5*^ in the effective model, with the spin 
scaled by a factor of one half. However, conservation 
of the total spin in the IHM does not lead to conserva- 
tion of total spin for the effective model, which is not 
S'C/(2)-invariant. In Table HIl we show the mapping of the 
most important quantities from the original ionic Hub- 
bard model to the effective spin-one model. 



C. Symmetries 

One relevant characteristic of the effective model is the 
extent to which the symmetries of the original model are 
preserved or modified. The interaction term is local, 
translationally invariant, and depends only on (S*^)^, in 
contrast to the on-site part of the IHM Hamiltonian in 
Eq. ([1]). The apparently greater translational symmetry 
of the effective model is a consequence of the reduction of 
state space in transforming to the effective model. Note 
that this is only true at half filling: the quantity (N) — L 
(see Eq. [7]) is conserved and breaks translational symme- 
try except at half filling, where it is zero. (Note that the 
interpretation of the ~ state is not translationally 
invariant.) Since the spin-exchange term has the same 



symmetries as the hopping term in the IHM, the remain- 
ing symmetries of the original model are preserved in 
the effective model. Conserved quantities in the origi- 
nal model, such as the total z-component of the spin, Sz , 
the total spin, s, and the number of particles, TV, are 
still conserved in the effective model, but have different 
meanings. 



II. NUMERICAL METHOD 

We have investigated the effective model by performing 
density-matrix renormalization group (DMRG) calcula- 
tions for different system sizes, from L = 200 up to 600 
sites, with open boundary conditions (OBC) and an even 
number of sitesi^>^ For small chains, boundary effects 
can be large, depending on the correlation length. Thus, 
in order to minimize any dispersion due to the edges 
Friedel oscillations and odd-even effectsj^^ we analyze 
systems of at least 200 sites. In order to achieve sufficient 
accuracy, at least 5 sweeps must be performed, with up 
to 1280 states retained in the last sweep. The maximum 
system size that can be accurately treated is then ap- 
proximately 600 sitcsi^ The maximum discarded weight 
of the density matrix for the effective model is always less 
than 10~*, and is typically zero to within the numerical 
precision far from the critical points In order to calcu- 
late ground-state properties, we target the ground state 
in the = sector; we target both the ground state 
and the first excited state in the = sector to calcu- 
late the 'cxciton' gap of the original IHM; and the lowest 
states in the S~ ~ 1 and 5*^ = 2 sectors are needed to 
calculate the charge and spin gaps, respectively, of the 
IHM 3 

We have repeated the same calculations using the 
dynamic block-state selection (DBSS) approach, fixing 
the threshold of maximum quantum information loss to 
X = 10~^ at each step i^^'^^ For instance, m w 500 basis 
states are enough to correctly describe the ground-state 
wave function of a system with 500 sites for e = 1.23. 
However, as we increase e the number of states required 
increases, for example to m 900 states for e = 2. For 
ground states of other symmetry sectors, e.g., the lowest 
triplet excitation, this number can sometimes be larger 
when the excited state is delocalized, despite the fact 
that its Fock subspace is smaller. Nevertheless, since 
we are interested in only the energy of these states and 
since measurements are carried out only on the absolute 
ground state, keeping of the order of a thousand states 
is usually sufficient. 

As the aim of the effective model is to describe the 
strong-coupling limit of the IHM when {U,A) >> t, we 
have compared results from the effective model to DMRG 
results for the IHM for ~ A = 20t^ All the quanti- 
ties that we measure: gaps, ionicity, bond order param- 
eter and polarization, are in agreement to within a few 
percent. 
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Table II: Mapping of relevant physical quantities to the effective spin one model. 



Quantity Ionic Hubbard Model Effective Spin One Model 



lonicity /= (^i,) ' = tJ2((^:] 



Polarization 

L ^ — i .tt ^ ' \ ^ " ^ / 2 

1 '"^^ 1 ^ 

Bond Order Parameter D = (-1)* (cjci+i + c+^jC;) D = y"(-l)' [ ( ( •S'i'' + 1 + -S.'^i + l) ■^i'+i 

-5= (s+sr^i + S-S++i))] 

L L 

AFM Order A = - 1)' (^r ) A=— ^(-l)'(sn 



III. 



BI TO SDI TRANSITION 



In this section we study the first transition between 
the band-insulator phase and the spontaneously dimer- 
ized phase. Wc have tuned the interaction coupling e 
starting from zero, where the system behaves like a band 
insulator, increasing it until the first transition point £ci 
is reached. In order to locate the transition point, we 
have studied the behavior of the singlet and triplet gaps 
and of the bond order parameter. The two gaps go to zero 
in the thermodynamic limit at the transition point and 
subsequently reopen. The value of the bond order param- 
eter, which measures the system's dimerization, changes 
from zero to a finite value across the transition. The 
existence of such a transition has been extensively dis- 
cussed for the muM^ Therefore, we have focused on 
the characterization of the transition by evaluating its 
critical exponents explicitly. 

The Hamiltonian of the effective model is new and 
is not evidently related to any known classical model. 
Therefore, we must first determine the value of the dy- 
namic critical exponent z in order to carry out finite-size 
scaling. Subsequently, we extract the correlation length 
exponent v from the divergence of the mass gap. Finally, 
the thermodynamic exponents (3, a and 7, which are all 
related to the free energy density, are obtained by ana- 
lyzing the divergence of the bond order parameter, the 
specific heat, and the bond-order susceptibility, respec- 
tively. 



of z can be different from one. Therefore, a determination 
of z is required to obtain and interpret all the remaining 
critical exponents. First, we identify the mass gap 



F{e,L)^E, {e,L)-Eois,L) 



(8) 



which is the gap that scales to zero most quickly close 
to the critical pointi^ This gap is proportional to 
where the correlation length ^ is limited by the system 
size L. Consequently, the ratio 

F{e,N) N 



R,{e,N,M) 



(9) 



F{e,M) M 

of the mass gaps for different system sizes behaves as 
Rz{ec,,N,M) - {N/M^'^ for N,M » 1, and thus 
depends only on the ratio of system sizes r = N/M^ In 



Fig. 2(a) we show that all the gap ratios with a particular 
r (r = 1.5 in the figure) cross each other at the same 
point, which is near Rz = 1. The behavior is similar 
for other values of r; we have examined r ~ 1.2, 1.25, 
1.33, and 2. In Fig. |2(b)[ one can see that curves with 
different r, scaled by the M = 200 gap, cross Rz = 1 
at the same point. Thus, it is clear that all curves cross 
each other at approximately the same value of e, e w 1.3, 
where Rz{N, M) w 1, consistent with z = 1j^ In order to 
carry out the scaling analysis of the critical coupling and 
other critical exponents, we take z = 1 in the following 
subsections. 



B. Correlation length exponent 1/ 



A. Dynamic critical exponent z 

For a quantum system related to a classical model by 
the transfer matrix, the dynamic critical exponent plays 
the role of an extra dimension, i.e., z = I. In general, 
space and time correlations can be coupled, and the value 



In order to proceed, we next need to calculate the criti- 
cal value of the coupling in the thermodynamic limit, Sci ■ 
The most efficient and accurate way of doing this is to 
carry out scaling using the logarithmic mass gap ratio,— 
defined as 

lnF{s,L + 2)-lnF{e,L) 



R{e,L) 



ln(L -f 2) - In L 



(10) 
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Figure 2: (Color online) Mass gap ratio Rz as function of the 
vario us system sizes A'' and M and the same 
(b) mass gap ratio for different r scaled by 



coupling e for (a) 
ratio r — 1.5, and 
the M = 200 gap 



This quantity can be used to define a sequence of pseudo- 
critical points for different system sizes using the criterion 
R{e*,L) + 1 = 0. In Fig. 3(a) wc show the behavior of 
i? (e*, L) + 1 as a function of e for various system sizes. 
The curves of the scaled ratio cross the line at two points, 
defining two sets of pseudo-critical points, which we des- 
ignate as £a{L) and el{L) for the lower and upper cross- 
ings, respectively. The finite-size scaling of bo th series of 
pseudo-critical points is depicted in Fig. 3(b)[ All curves 
are fit with third-order polynomials in 1/L. In the ther- 
modynamic limit, e* and el converge to the same point 
to within the accuracy of the extrapolation, confirming 
that the transition is second order. The finite-size scaling 
of the position of the minimum in the mass gap provides 
an alternate way of determining 6^ ■ This can either be 
done using the mass gap, Eq. ([5]), directly, which we des- 
ignate as e,„(L), or using the minimum of the mass-gap 
ratio, Eq. 0, designated as £r{L)- The extrapolations of 
positions of the minima, Em and £r also converge to the 
same point, providing a confirmation of the consistency 
and stability of the extrapolation procedure. We obtain 
the location of the critical point at 



£c, = 1.286(5). 



(11) 



We can now estimate the correlation-length exponent 
using the finite-size version of the Callan-Symanzik (3- 
function^i^^^ 



1 dF{£,L) 



(12) 



F (e, L) de 
which has critical behavior 

fi,,{s,,,L)^ L-^ . (13) 
To calculate the exponent v, we proceed as follows: Given 



Figure 3: (Color online) (a) Logarithmic mass gap ratio as 
a function of e for various system sizes and (b) finite-size 
extrapolations of the critical point using both sequences of 
pseudo-critical points, as well as the two definitions of the gap 
minimum. Here e,„ is the position of the mass gap minimum, 
Et the position of the mass-gap ratio minimum, and e'a and el 
are the upper and lower sequences of pseudo-critical points, 
respectively. The lines are guides to the eye. 



a sequence of pseudo-critical points, e*{L), we extrapo- 
late the ratio of the /3-functions for different system sizes 



/?cs {e*,L + i 
Pes {e*,L) 



L + l 
L 



to the thermodynamic limit. Here it is important to 
choose e*{L) carefully: extrapolating using a series of 
pseudo-critical points that is close to the gap minimum 
can yield unreliable results because the derivative of the 
mass gap remains zero or close to zero. Therefore, we 
utilize the ratio from the second scries of pseudo-critical 
poi nts = rather than from the first e'^{L) [see 

Fig. [3(b)]. If ll/v] < 1 and L » i, then 



L 

1 



- 1 



Pcs{e*,L) 

From the numerical extrapolation, we obtain 
- = 0.996(5) . 



(14) 



A plot of unsealed mass-gap data and its collapse using 
this scaling exponent is shown in Fig. [4] 



C. Thermodynamic exponents P, a, 7 

The bond order parameter characterizes the bond- 
order-wave (BOW) phase. Fabrizio, Gogolin, and Ners- 
esyan have argued that the bond order parameter is the 
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Figure 4: (Color online) Scaling of the mass gap around the 
first crit ical point: (a) unsealed data. The lines are guides to 
the eye. |(b)| Rescaled data F (e, L) L are plotted as a function 
of the rescaled coupling L (e ~ Eci ) /eci . 



Figure 5: (Color online) Bond order parameter as a function 
of the coupling e around the transition point: (a)| data for 



different system sizes near the first transition and (b) data 
rescaled as D (e, L) L^^^ plotted as a function of the rescaled 
coupling L (e - 6^ ) /sci . 



right quantity to characterize the Ising transition in the 
IHM.HiiS The order parameter, expressed in the spin-one 
language, is given by 



1 



L - 1 ^ 

i=l 



Si I S: S, 



Sa s. 



l+l 



(15) 



DMRG results for the bond order parameter as a func- 
tion of the coupling e near the first transition point are 



depicted for various system sizes in Fig. 5(a) 



We can use the bond order parameter to determine the 
associated critical exponent /3, i.e.. 



D{£ ^ £ci,L) ~ L^^ . 

Using the logarithmic derivative 

\nD {e\L + t) - InD {e\L) 



we obtain 



In {L + ()- In L 



- = 0.124(5). 

V 



1 

V 



(16) 



(17) 



(18) 



The excellent data collapse of the rescaled data, as can 



be seen in Fig. 5(b) confirms that the transition point 
belongs to the 2D Ising universality class. Results for the 
finite-size scaling of the exponent (3 are plotted in Fig. [6l 

Since, in a quantum phase transition, the coupling 
plays the same role as temperature in a thermal phase 



0.5 - 




-0.5 - 



0.001 0.002 0.003 0.004 0.005 

1/L 

Figure 6: (Color online) Finite-size behavior of the exponents 
V, (3, a, and 7. The fit is to a third-degree polynomial in 1/L. 
Note that the points for the two smallest 1/L are not included 
in fitting 7. 



transition, we can define a corresponding "specific 
heat"^M2 



C.„ (£,i) = -- 



T 



Note that this quantity docs not correspond to the real 
specific heat. Nevertheless, due to the scaling relations 
and its interplay with the other quantities, it has to di- 
verge with the exponent a. The physical specific heat 
exponent is related to our a by the Griineisen parame- 



ter 



43 



The specific heat usually contains a regular term that 
is typically larger in amplitude than the singular one. 
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Therefore, instead of using the fogarithmic derivative to 
estimate the exponent ajv, we instead use the ratio 



L Ct, (e, i + 2) — c„ (e, V) a 
2 Cv (e, L) V 



(19) 



To overcome possible problems in determining this expo- 
nent, we use the Hcllman-Feynman^ theorem to exploit 
the accuracy of the DMRG in calculating local quantities 



de 



1 



(20) 



This trick reduces the computational cost to that of cal- 
culating the first derivative of the cubic spline, which 
interpolates the data points 4^ The result is the follow- 
ing: 



a 



0.00(1). 



(21) 



The finite-size behavior of the various exponents is plot- 
ted in Fig.[6l The scaling relation a = 2 (1 — i/) is fulfilled 
by Eqs. ^ and 

Finally, we determine the exponent 7 associated with 
the relevant susceptibility. The susceptibility corre- 
sponding to the bond order parameter is 



XD (e,i) = - 



1 dD{e,L) 
L 



dh 



D 



(22) 



In order to calculate this quantity, we turn once more 
to the Hellman-Feynman theorem and to linear response 
theory. We perturb the Hamiltonian with a small field 
hn conjugate to the order parameter D. The field has to 
be small enough to reveal a linear regime in the changes, 
but not smaller than the actual DMRG resolution; we use 
2dhD = W~H. We have measured the order parameter 
for four points around hu = in order to compute its 
first derivative at h]j = 0. 

Once we have evaluated the static susceptibility for 
different system sizes, we proceed in the same way as for 
the previous exponents. The scaling relation is 



Thus, from 



IiiXd{£*,L + £) -\nxDie*,L) 7 



In {L 



InL 



(23) 



(24) 



we obtain the last thermodynamic exponent, as plotted 
in Fig. [6l with the value 



- = 1.72(5). 



(25) 



As shown in the figure, the last points for the largest 
system sizes have been excluded in calculating the expo- 
nent. The reason is that the calculation of the suscep- 
tibility becomes uncontrolled for very big system sizes. 
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L ( t- ) / 
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Figure 7: (Color online) (a) Bond order parameter suscepti- 
bility as function of the coupling e for different system sizes, 
(b) The collapsed curves scaled using the exponent 7 = 7/4. 



In order to compensate the occurrence of nonlinear be- 
havior in the response for larger system sizes, we would 
have to use a very small perturbation field. However, the 
effect of such a small field can be difficult to distinguish 
from the numerical noise. In addition, we have to carry 
out two cubic-spline interpolations: one to determine the 
derivative of the bond order parameter as function of the 
perturbation field and one to fit its susceptibility. For 
these reasons we neglect the points at the two largest 
system sizes. We see that the second scaling relation 
^ = 2{v — (3) is fulfilled to within our estimated error i22 
Other quantities, such as the electric polarization and the 
electric susceptibility, scale with the same exponents as 
the bond order parameter and the bond-order suscepti- 
bility, rcspcctivelyf^i 

In addition, we have calculated the value of the central 
charge governing the underlying conformal field theory 
numerically in two different ways. In the first method, 
we use that the scaling of the low-lying energy levels 
with system size is uniquely determined by the confor- 
mal tower4^ This scaling can be used to determine the 
central charge4i The value obtained, c = 0.50(4), is con- 
sistent with that expected for the 2D Ising model. In the 
second method, we determine the central charge from 
the entropy profile, which has a known form dependent 
only on the central charge4^ We obtain the same value 
(c 0.5) to within the numerical accuracy at 8^- 



IV. 



SDI TO MI TRANSITION 



In this section, we present numerical results on the sec- 
ond transition where the system passes from the sponta- 
neously dimerized phase to the Mott insulator phase with 
increasing e. We will show both how the spin gap closes 
when approaching the critical point from below and 
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how the bond susceptibiUty diverges when approaching 
£c2 from above. Our results confirm the KT scenario 
with an essential singularity at e^^ and a critical phase 
for e > ec2- In order to do this, a more careful treatment 
than at the 2D Ising transition point is required. 



A. Correlation length and mass gap 

For a KT transition, the correlation length diverges 
exponentially as the transition point is approached from 
the gapped phase and remains infinite in the critical re- 
gion that follows 4^ Since the mass gap is related to the 
inverse of the correlation length, the mass gap has to 
close exponentially as the transition is approached and 
is zero in the critical region. However, for finite-size sys- 
tems, the correlation length ^ is limited by the system 
size L. Very large system sizes or the inclusion of higher 
order corrections are required to reveal the exponential 
divergence, which is restricted to a narrow region close 
to the KT transition. 

In order to locate the position of the second transition 
point £c2j we analyze the scaling of the mass gap, de- 



picted in Fig. 8(a) The finite-size scaling analysis made 



for the first transition cannot be used here because suf- 
ficiently large systems to study the logarithmic scaling 
cannot be reached. Instead, we prefer to use a different 
approach based on conformal field theory (CFT). Within 
the Mott insulator phase, where the spin sector is gap- 
less, the system is critical and can be described by a CFT. 
Furthermore, the characteristic excitation gaps scale with 
system size L as 



E, (L) - So [L] = 



2TTXiV 



(26) 



where Xi is the corresponding scaling index and v is the 
"excitation" velocity. Since this expression is valid only 
in the critical region corresponding to the Mott insulator, 
the extent to which it is fulfilled can be used to locate 
the transition point. In a plot of the mass gap times 
the system size L, Fig. 8(b) all curves merge into a sin- 



gle one exactly at a critical point , as expected from 
Eq. (pS)) . Therefore, the system is in a critical regime 
above a critical coupling 



1.8(1), 



(27) 



The point at which the curves merge is clearly separated, 
see Eq. ([Tl]) . from the first critical point that we found 
at • 

An analysis of the mass gap ratio, sec Eq. PT!]) . is also 
useful. In contrast to what happens at the first transi- 
tion, , the curves do not cross the line corresponding 
to a ratio of unity due to the logarithmic corrections4i 
Nevertheless, the curves remain very close to zero ev- 
erywhere in the critical region above , as can be seen 



in Fig. 9(a) In the region preceding Ccg, the value of 
the mass-gap ratio increases with the system size, as ex- 
pected for a gapped system. The overall behavior of the 
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Figure 8: (Color online) (a) Mass gap and (b) mass gap times 
L relative to the second critical point as function of the cou- 
pling e and for different system sizes. 



mass-gap ratio curves further confirms that there is a sec- 
ond transition point at and supports the KT scenario. 

In addition, we define and calculate the scaled differ- 
ence of mass gaps, Q, 



Q{s;L',L) = 



U F{e,U) ■ U - F{e,L) ■ L 



2n 



U -L 



(28) 



For an arbitrary V, the first-order finite-size scaling 
terms cancel out and Q{e;L',L) vanishes in the criti- 
cal region. In Fig. |9(b)| we show results for L' = 500. 
We conclude that the second critical point occurs at 
~ 1-8 and the gap closes exponentially. Since the 
distance between the two critical points is much bigger, 
£c2 ~ £ci ~ 0.5, than any deviation due to the logarith- 
mic corrections, we conclude that there are two phase 
transitions. 

We have also calculated the approximate /3-function 
Pes, Eq. However, for this kind of transition, it 

has no zeros (as expected) 4i Nevertheless, we can ex- 
trapolate the value of the minima of the /3-function as a 
function of the system size to the thermodynamic limit. 
This yields an alternate estimate of ec2j =1.9(1). 



B. The bond-order and electric susceptibility 

In order to classify the transition as a KT transition, 
we examine the bond-order susceptibility and the electric 
susceptibility, see Fig. [TOl The behavior of the peak of 
the bond-order susceptibility can be used to estimate the 
exponent of the susceptibility. 



7peak {L) 



peak 



L + 2 



peak' 



In (L-h2) -InL 



(29) 
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scenario of an infinite-order phase transition. The criti- 
cal exponent of the susceptibiUty 7 cannot be determined 
accurately because of strong finite-size effects due to the 
strong influence of the bond-order wave which scales to 
zero very slowly, i.e., as 

Additionally, we have made a preliminary calculation 
of the central charge from the entropy profile assuming 
that it has the form predicted by CFTi^ In order to de- 
termine the transition point, we minimize the ^(^ of the 
fit to the conformal form and confirm that c w 1 (c = 1 is 
expected for this type of KT transition) at this point. We 
obtain a rough estimate of £c2: ~ 1.65(15), which is 
consistent with the results of our finite-size scaling anal- 
ysis, = 1-8(1) to within the accuracy of the scaling. 
Thus, the three estimates of the critical coupling, £c2i 
e^^ , and e^^ are consistent with one another. Our best 
estimate is given by , since the other two estimates are 
rougher and more likely to contain systematic errors. 



(b) 



Figure 9: (Color online) (a) The logarithmic mass gap ratio 
plus unity, (b) The scaled difference of mass gaps Q for L' — 
500 and various values of L. 
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Figure 10; (Color online) (a) The bond-order susceptibility 



and (b) the electrical susceptibility for the SDI-MI transition. 



The position of the peak in the bond-order susceptibil- 
ity converges to the value epoak ~ 1.62, and the series of 
pseudo-exponents, 7peak (^)' converges to 7 w 1.27 in the 
thermodynamic limit. For comparison, we calculate the 
electric susceptibility, shown in Fig. |10(b)[ The finite- 
size effects are much stronger for the electric suscepti- 
bility than for the bond-order susceptibility. In fact, we 
also observe a narrow peak in Xe that grows and moves 
with the system size. In general, we conclude that the 
coincidence of the mass gap closing to zero exponentially 
and a diverging susceptibility corresponds to the typical 



V. DISCUSSION 

We have analyzed the band-insulator-to-Mott- 
insulator transition in the strong-coupling limit. Using 
simple strong-coupling arguments, we have derived an 
effective model starting from the ionic Hubbard model. 
The effective model, which we have formulated in a 
spin-one representation, captures the physics of the 
transition and is less computationally demanding than 
the ionic Hubbard model. It contains spin-exchange 
processes which are strongly restricted compared to 
those of a conventional spin model. In addition, the 
effective model demonstrates that a single interaction 
parameter governs the transition. Our density-matrix 
renormalization group study of this model confirms 
that there are two transitions at two clearly separated 
coupling strengths. The system undergoes a transition 
from a band insulator to a spontaneously dimerized 
insulator followed by a transition from the spontaneously 
dimerized phase to a Mott insulator with increasing 
effective interaction. This behavior corresponds to the 
behavior of the ionic Hubbard model found in previous 
work. 

In Fig. 1111 we explicitly compare the phase boundaries 
obtained in our work to phase boundaries obtained nu- 
merically for the ionic Hubbard model in Refs.[§ andfiol. 
The phase diagram is plotted in the 45° rotated U-A 
plane of the ionic Hubbard model, so that the abscissa 
corresponds to our effective parameter e = U — A and 
the intermediate phase is expanded relative to the de- 
piction in Fig. [1] Since our effective model is based on a 
strong coupling expansion in U and A, our results should 
be applicable to the ionic Hubbard model in the large 
U + A limit. As can be seen, for the BI-SDI bound- 
ary, both ionic Hubbard model results tend towards our 
strong-coupling value as U + A becomes larger, although 
the largest coupling point (at A = 20) from Rcf. 10 is 
still outside our error bars. The results for the SDI-MI 
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transition boundary have larger discrepancies, but our 
results lie between strong U + A extrapolations of the 
phase boundaries of Ref. [13 and that of Ref. [^. This 
underlines the difficulty of obtaining the transition point 
in this infinite-order Kosterlitz-Thouless transition. To 
within large, but realistic error bars, the three sets of 
results for this phase boundary arc not necessarily in- 
consistent with each other. 
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Figure 11: (Color online) Rotated ground-state phase dia- 
gram phase diagram of the ionic Hubbard model depicting 
the phase boundaried obtained in Refs. |^ and [l^, as well as 
the transition points Sci = 1.286(5) and = 1-8(1) obtained 
in this work, which apply in strong coupling in U /t, A/t. The 
estimated error in our results are indicated by the gray-shaded 
bars at the upper and lower axes. 



excited state, is the mass gap, and it goes to zero at , 
while the spin gap (the gap to spin triplet excitations) 
remains finite. For < s < mass gap is set 

by the gap to the lowest-lying triplet, i.e., the spin gap, 
which is degenerate with singlet excited states. This gap 
goes to zero at ■ For £ > > the spin and cxciton gaps 
remain zero, as expected in a critical phase, but the gaps 
to add or remove one or more particles remain finite. 

We note also that the mapping of electronic systems 
to spin-one systems derived here can be adapted to a 
larger class of similar models or to generalizations of 
the ionic Hubbard model, e.g., to chains with an ionic 
potential with a different periodicity or even to two- 
dimensional systems. Another potentially interesting ap- 
plication would be to relate exactly solvable spin-one 
models to electronic models and vice versa via the spin- 
one composite representation!^ 
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Our extraction of the critical exponents for the first 
transition confirms that it belongs to the two-dimensional 
Ising universality class. We have also shown that the uni- 
versal scaling relations are fulfilled to within our numeri- 
cal accuracy. At the second transition, we have observed 
that the mass gap closes exponentially and that all rel- 
evant susceptibilities diverge. An analysis of the scal- 
ing of the mass gap and of the bond-order susceptibility 
confirms typical Kosterlitz-Thouless behavior. The mass 
gap closes at the critical point and then remains zero 
as the interaction is further increased. The susceptibil- 
ity diverges in the entire critical region above the second 
transition point. The overall scenario, with an Ising-like 
transition from the band insulator to the spontaneously 
dimerized insulator followed by an infinite-order transi- 
tion from the dimerized insulator the Mott insulator, is in 
complete agreement with the field-theoretical prediction 
for the ionic Hubbard modeliii 

The evolution of the appropriately mapped gaps with 
increasing e in the effective model is consistent with the 
picture obtained for the ionic Hubbard model in Ref. [Tol . 
Deep in the band insulating phase, all gaps in the spin 
and charge sectors are equal and are set by the band 
gap. As is approached, the exciton gap, defined as the 
energy gap between the ground state and the first singlet 



Appendix: DERIVATION OF THE EFFECTIVE 
MODEL 

The effective Hamiltonian can most easily be derived 
by first expressing the original Hamiltonian as a func- 
tion of the Hubbard operators X"^ = \ai){(3i\, where the 
\ai) and designate an element of the Hubbard ba- 
sis {|0), I t), I i), \d)} on site i. Exphcitly, they can be 
expressed as 

"(i - np(i - fit) ct(i-nj.) Ci(i-nT) cic^ ' 

_ cj(i-nx) (i-fix)ni c|c| -cih^ 

^ c|(i-n|) cjc| n|(i-n|) n|C| 

c|c| ~c|n| fiicj fiifi-^ 

For instance, we rewrite the ionic potential and the 
Coulomb interaction as 

L 

Hu = C/^Xf 

i=l 
L/2 
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and 



becomes 



Ha 



X. 



TT 



X. 



ii 



dd 



2 ^ 

i=i 



2X, 



2X 



They can then be mapped onto the spin-one model ex- 
pressed in terms of the operators if* = |si)(s-|, with 
\si) the spin-one Sz basis {|1), |0), | — 1)} on site i. The 
single-site Hilbcrt space truncation is defined as 



X. 
X. 

X. 



Q/3 

i 

a/3 
i 



> for a OT (3 ~ and i 
y for a or f3 — d and i 
L"^ otherwise. 



2j 
2j 



(A.l) 



In the spin-one basis, 
'{sff+s} 



L = 



2 

S~ SI 



V2 



V2 



(A.2) 



Hence, the interaction and the potential parts are trans- 
formed to 



Ha 



Hn 



-T 

2 ^ 

-t 



L/2 



00 

2j-l 



(A.3) 



L 



-1-1 
2i-i 



2L 



2j-l 



2j 



L. 



2j 



(A.4) 



Altogether, defining the coupling constant e — C/ — A, 
the doping 5 = N — L, and writing the terms using spin- 
one operators, see Eq. ()A.2|) . the two-term contribution 



i=l 



2 2 



(A.5) 



Likewise, the hopping part is translated to 

H^ = tj2(Lr'i^U-LTL%l 

1=1 

or, in the spin one language 



(A.6) 



i=l 



(A.7) 



which is equivalent to Eq. (O. A sketch of the allowed 
processes is shown in Fig. [T^l These processes arc a rela- 
tive small subset of those of the isotropic Hciscnbcrg spin 
chain model. Note that the AFM exchange in the IHM 
maps to a sequence consisting of two scattering processes 
in the effective model. 




Figure 12: (Color online) Sketch of the allowed processes, 
which are a relatively small subset of those of the isotropic 
Heisenberg spin-chain model. 
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